/* Author: Panos Sakkos 2011 
 * Email: p.sakkos@di.uoa.gr
 */
 
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <iomanip>
#include "trigonal.h"
#include "vector.h"

#define THREE_DIGITS_PRECISION 0
#define PRINT_X_VECTOR 0

using namespace std;

int main (int argc, char *argv[])
{
	if(argc != 4)
	{
		cerr << "Wrong arguments!\n";
		return EXIT_FAILURE;	
	}
	
	/* get n, desired accuracy of solution and max iterations threshold from the parameters */
	
	int n = atoi(argv[1]);
	double accuracy = 0.5 * pow(10.0, (-1) * atoi(argv[2]));
	int iterationThreshold = atoi(argv[3]);	
	
	Trigonal trigonal(n);	

	/* Read non-zero input for Trigonal object */

	trigonal.ReadInput();
	
	Vector b(n);
	cout << "Please, type input for b vector: \n";
	b.ReadInput();
	
	Vector x_previous(b);
	Vector x_next(n);
	
	/* For t = 0.1, 0.2, ..., 1.9 */

#if THREE_DIGITS_PRECISION
	for(float t = 0.100; t < 2.000; t += 0.001)
#else
	for(float t = 0.1; t < 2.0; t += 0.1)
#endif
	{
		int iterationCount = 0;
		bool approximationFound = false;
	
		while(iterationCount < iterationThreshold)
		{
			/* For every value of x_next vector */
			
			for(int i = n - 1; i >= 0; i--)
			{
				double USumNext = 0.0;
				double USumPrevious = 0.0;
				double LSum = 0.0;
								
				for(int j = i + 1; j < n; j++)
				{
					USumNext += (trigonal.GetValue(i, j) / trigonal.GetValue(i, i)) * x_next.GetValue(j);
				}	
								
				for(int j = i + 1; j < n; j++)
				{
					USumPrevious += (trigonal.GetValue(i, j) / trigonal.GetValue(i, i)) * x_previous.GetValue(j);
				}

				for(int j = 0; j < i; j++)
				{
					LSum += (trigonal.GetValue(i, j) / trigonal.GetValue(i, i)) * x_previous.GetValue(j);
				}
				
				double value = (1-t) * x_previous.GetValue(i) - USumNext - (t - 1) * USumPrevious - t * LSum + t * (b.GetValue(i) / trigonal.GetValue(i, i));
				x_next.SetValue(i, value);
			}
			
			float norm = (x_next - x_previous).GetInfinityNorm();

			if(norm < accuracy)
			{
				cout << "Vector approximation for t = " << setfill (' ') << setw (7) << t << " was found ";
				#if PRINT_X_VECTOR
				x_next.Print();
				#endif
				cout << " after " << iterationCount + 1 << " iterations\n";
				approximationFound = true;
				break;
			}
			else
			{
				x_previous = x_next;	
				iterationCount++;
			}
		}
		
		if(approximationFound == false)
		{
			cout << "No approximation found for t = " << t << " in " << iterationCount << " iterations" << endl;
		}
		
		approximationFound = false;			
	}
		
	return EXIT_SUCCESS;
}

